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Abstract 

We study the dilute Fermi gas at unitarity using molecular dynamics with an effective quantum 
potential constructed to reproduce the quantum two-body density matrix at unitarity. Results for 
the equation of state, the pair correlation function and the shear viscosity are presented. These 
quantities are well understood in the dilute, high temperature, limit. Using molecular dynamics we 
determine higher order corrections in the diluteness parameter raA 3 , where n is the density and A 
is the thermal de Broglie wave length. In the case of the contact density, which parameterizes the 
short distance behavior of the correlation function, we find that the results of molecular dynamics 
interpolates between the truncated second and third order virial expansion, and are in excellent 
agreement with existing T-matrix calculations. For the shear viscosity we reproduce the expected 
scaling behavior at high temperature, rj ~ 1/A 3 , and we determine the leading density dependent 
correction to this result. 
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I. INTRODUCTION 



There has been a significant amount of interest in equilibrium and transport properties 
of cold atomic quantum gases near a Feshbach resonance MiM- The Feshbach resonance 
is used to tune the energy of a molecular bound state to the threshold of a scattering 
state, causing the scattering length to diverge. At resonance the scattering cross-section is 
universal, independent of the microscopic details, and is limited only by unitarity. In this 
regime there are no small parameters that can be used to justify a perturbative expansion. 
The only ab-initio approach to the problem is the quantum Monte-Carlo (QMC) method 



llj . QMC simulations have been used successfully to compute thermodynamic properties, 



but it is very difficult to use imaginary time QMC simulations to compute real-time response 
functions and transport coefficients (see, however, [3] for a recent attempt). This implies 
that we lack reliable bench mark calculations for experimental attempts to determine the 
viscosity of a degenerate Fermi gas near unitarity. While the viscosity can be computed 



reliably at both high |l3Hl5j and low [if], [l^ temperature there is no systematic approach 



in the strongly coupled regime T ~ Tp, where Tp is the Fermi temperature. A number 
of authors have used diagrammatic methods to study transport properties in this regime 
but it is not a priori clear what kind of diagrams have to be included. 
In this work we introduce a novel approach to the dynamics of quantum gases. The 
method is based on a classical molecular dynamics simulation in which quantum effects are 
encoded in an effective classical interaction among the atoms. The quasi-classical iV-body 
interaction is constructed such that the classical calculation exactly reproduces the diagonal 
component of the quantum iV-body density matrix. This implies, in particular, that the 
iV'th quantum virial coefficient is reproduced. In this paper we restrict ourselves to two-body 
terms in the interaction. In this case the second virial coefficient will be reproduced exactly, 
and through molecular dynamic simulations, the one and two-body components of the higher 
virial coefficients are resumed. The strength of the molecular dynamics method is that very 
complicated many body correlations are taken into account. The drawback is that genuine 
quantum many-body effects such as pairing and superfluidity are not included. Quasi- 
classical molecular dynamics has been used successfully in the study of strongly correlated 



Coulomb plasmas 
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211 ] . However, our work differs in that the quasi-classical potential 



for unitary fermions is of purely quantum mechanical origin. In contrast, the quasi-classical 



2 



Coulomb potential, known as the Kelbg potential 22j, is a small correction to the classical 
1/r behavior. 

This paper is organized as follows. In Sect. HT1 and IHII we introduce the method and derive 
the quasi-classical potential at unitarity. In Sect. IIVI we describe the molecular dynamics 
simulations and in Sect. |V]we present results for the equation of state, the pair correlation 
function, and the shear viscosity. 

II. THE PARTITION FUNCTION 

In this section we introduce a classical partition function which is equivalent, order-by- 
order in a cluster expansion, to the full quantum partition function. The classical partition 
function depends on /c-body potentials which can be determined from the quantum mechan- 
ical Slater sums. In this work we will restrict ourselves to the case k = 2, which means that 
the second virial coefficient is reproduced exactly. We will follow the notation used in 

The quantum mechanical partition function for a system of N particles is given by 



Zn ~ N\X" 



^ J {dv l ...dv N )W^{v 1 ,...,v N ) , (1) 
where we have defined the iV-particle Slater sums 

^ (r 1 ,...,r N ) = N\X 3N |*« (ri, ...,**) . (2) 



In this expression A = y mk^T * s ^ ne thermal wavelength, j3 — 1/T is the inverse tempera- 
ture, and ty a ( r i, ■ ■ ■ , r Jv) is the wave function of an N particle state with energy E a . The 
partition function can be expanded systematically in powers of n\ 3 in terms of the virial 
coefficients bi, 

where {mi} is a set of integers mi > that satisfies the constraint YliLi^ m i = N. The 
corresponding expression for the pressure has the form 

oo 



p 1 



kT X 3 

i=i 
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where z = is the fugacity and /i is the chemical potential. Each virial coefficient can 
be expressed in terms of integrals over the iV-particle Slater sums. For example 

b 2 = ^7 J dv 1 dv 2 [^ (2) ( ri ,r 2 ) - W {1) (r^W^ (r 2 )] , (5) 

6s = 3!W / dridT2dr3 ^ (3) (ri,r 2 ,r 3 )-^(r 1 ,r 2 )^ 1 )(r 3 ) 

- W^(r 2 ,r 3 )W^\r 1 )~W^(r 3 ,r 1 )W il \r 2 ) + 2W {1 \r 1 )W^\r 2 )W il \r 3 )] . (6) 

We may compare these results to the corresponding expressions for a classical system. We 
consider the most general partition function containing arbitrary iV-body interactions 

Z N = J {dr 1 ,...,dr N )e-^<i^-^<i<^" + - . (7) 

The virial expansion of the classical partition function has the same form as the quantum 
expansion in Eq. fl3]), but the expressions for the virial coefficients are different. We have 

63 = —L- J dndr 2 dr 3 [ e -^^+vi2+v 23 +v 13 ) _ e -pv 12 _ e -0v 23 _ e -Pv vi + 2 ] _ ( 9 ) 

It is clear that one can construct a classical iV-body potential so that the classical and 
quantum virial coefficients agree order by order. For example, we can construct effective 2 
and 3-body potentials 



-p-Hogtyv® (Ti,^)) , (10) 

_i W (3) (ri,r 2 ,r 3 ) 
v m = -P log w(2) j- w(2) j- - ^ w(2) j- - - j , (11) 

such that the classical system described by the partition function in Eq. (J7j) will have the 
same virial coefficients as the quantum mechanical system at the same order. In practice 
we will truncate the expansion at second order. We note that even in this case we retain 
all contributions to the third and higher virial coefficients that arise from powers of W®. 
Only genuine three and higher-body correlations are missing. 



III. QUASI-CLASSICAL TWO-BODY POTENTIAL AT UNITARITY 



As discussed in the previous section the first non-trivial virial coefficient can be repro- 
duced by introducing the effective two-body potential defined in Eq. (FID]). The potential 
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depends on the logarithm of the two-particle Slater sum W^ 2 \ This quantity is defined in 
terms of the solutions of the two-particle Schrodinger equation 

mi* (ri,r 2 ) = E a ^ a (n,r 2 ) , 

W = -^(Vi + V^)+V(|r 1 -r 2 |) . (12) 

In the case of a spherically symmetric interaction the two particle Slater sum is only a 
function of the relative coordinate r = |r 2 — ri|. We find 

W (2) (r) = 2 5/2 A 3 y (2?_M) y R kl {rfe-^ , (13) 

Z — / 47]- Z — <■ 

where Rki(r) satisfies the radial Schrodinger equation 



i?w(r) = 0, k 2 = ^6, . (14) 



A. Free particle 



It is instructive to begin by deriving the effective classical potential for a free particle. 
This potential takes into account the effects of quantum statistics: repulsion for identical 
fermions, or attraction for identical bosons. The wave function of a free particle is 

R kl (r) = V2^ 3l (kr) , (15) 

and the corresponding two-particle Slater sum is 

(r) = 2 7 / 2 A 3 J2 / k " dk Ji(kr) 2 e-^ . (16) 

We compute the Slater sum for identical bosons and fermions by restricting the sum over 
all states to I = even or I = odd, respectively. We find 

ol/2\3 r 22 2 2 

W {2) (r) = — — / dk k 2 [1 ± sinc(2A;r)] e~~ = i ± e — fr- ? ( i7 ) 

where sinc(:r) = sm(x)/x and we have made use of the identities 

E/^, ,\ -2/, \ 1 + sinc(2A;r) 
(21 + l^ftAx) = ^ ^ , 

Z=0,2,4,- 

E. 9 . 1 — sinc(2/cr) . . 

(21 + l)jf (At) = f . (18) 

i=l,3,5,... 

The resulting potentials are therefore 

u^f = -k B T log (l ± e - 2 - 2 /A 2 j ; (19) 

with the '+' sign for bosons and the '-' sign for fermions. 



B. Classical potential at unitarity 



At unitarity the physics is independent of the precise form of the interaction potential. 
We will therefore treat the interaction among opposite spin particles as arising from an 
attractive square well of depth Vq and range b. In the zero range limit only s-wave scattering 
contributes to the scattering amplitude. Outside the range of the potential the I = wave 
function has the form 

R t ,^ = VW Sia{k l +S ° ) , (20) 

where 5q is the s-wave phase shift. The correction to the free I = Slater sum is given by 

AWffi = [ dk [sin 2 (At + 5 (k)) - sin 2 (£r)l e~^r . (21) 

V2r 2 J 

In principle analogous expressions can be found inside the interaction region r < b. However, 
at unitarity, we are interested in the limit b — > while keeping y/V^b = n/2 fixed. We can 
therefore evaluate Eq. (J2"T]) for 5 = it/ 2 independent of k. In this case the integral is 
straightforward. We find 

= ^e—*A* , (22) 



and the effective potential at unitarity is given by 



A 2 



<T = -k B T\og 1 + — e" 2 ^ 2 . (23) 



We observe that the only length scale in the potential is the thermal wave length. 



IV. MOLECULAR DYNAMICS SIMULATION 



Having constructed the quasi-classical effective potential we now describe the molecular 
dynamics simulations. We consider a two component system with N = Nf + N± particles. 
The net polarization is zero and Nf = jVj, = N/2. The two body interactions between like 
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and unlike spins are 1 

u n = u;t = ~k B T log (l + ^- e -^iy? 
un = u u = -k B T log (l - e - 2 - 2 /A 2 j (24) 

The molecular dynamics equations of motion are 

f = ^, f = F, (25) 
at m at 

where % = 1, . . . , N and Fi is the force on the i'th particle due to the potential given in 
Eq. We measure the temperature in the simulation from the average kinetic energy 

per particle, 

r =s?(s> (**)')• < 26 > 

where the angular brackets denote an average over the simulation time. The pressure is 
computed using the virial theorem 



PV = nT+~(^r i -F-\ 



(27) 



It is advantageous to adopt a system of dimensionless units in which to perform the molecular 
dynamics simulations. We have used the system of units described in Table HI In particular, 
we use the thermal wave length A as the unit of distance, and A(m/T) 1//2 as the unit of 
time. We will denote quantities that are expressed in simulation units by a star, for example 
r* = r/X. In this system of units the simulation temperature T* is equal to unity, and the 
physical temperature is adjusted by changing the density n* = nX 3 . The ratio T/Tp, where 
Tp = (3Tc 2 n) 2 / 3 h 2 1 (2mkB) is the Fermi temperature, is given by T/Tp = 47r(37r 2 n*)~ 2//3 . 

We note that the effective quasi-classical potential is temperature dependent. In practice 
we choose a simulation density n* = nX 3 . We initialize the simulation using a guess for 
the total kinetic energy. We monitor the kinetic energy as the system equilibrates and add 



1 In practice one must soften the singular behavior of the attractive interaction between opposite spins. For 
this purpose we have replaced the term 1/r 2 inside the logarithm in Eq. (|24|) by (1 + \/~2nlo) / (r 2 + Iq), 
where is a regularization parameter. The normalization was chosen so that the second virial coefficient 
62 is insensitive to changes in Iq- We have checked that our numerical results are insensitive to Iq within 
the errors quoted in the text as long as Iq < 0.05 in the simulation units defined below. 
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Simulation Units 



Mass 


m= 


=mass of one atom 


Length 


A 




Energy 


k B 


x system temperature 


Time 


t* -- 


= Xyjm/T 


Derived Units 


density 


n* 


= NX 3 /V 


Temperature 




= 1 


Pressure 




= PX 3 /T 


Shear viscosity 


rj* 


= r/A 3 /Tt* 



TABLE I: System of units used in the numerical simulations of this work. 

or subtract energy by rescaling the velocities to reach the desired simulation temperature 
T* = 1. After the system equilibrates at this temperature we measure observables like the 
equation of state, the pair correlation function, and the correlation function of the stress 
tensor. The main simulation parameters are listed in Table [III Fluctuations are used to 
determine statistical errors. In addition, there are a number of systematic errors whose effect 
are difficult to quantify. One such systematic error is due to the finite number of particles. 
For the equation of state we have performed calculations with A" = 32, A" = 108 and A" = 256 
particles. Fig. [1] shows that the corresponding results agree within the statistical errors. Due 
to computational limitations we have only used runs with A" = 108 for our calculations of 
the pair and stress correlation functions. Another source of systematic uncertainty is the 
fact that molecular dynamics simulations are most naturally performed in a microcanonical 
ensemble (at fixed energy). This requires us to tune the energy very precisely in order to 
reach the simulation temperature. This is challenging because of long equilibration times 
and finite size fluctuations. 
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Simulation Parameters 

^atoms 108 

15 0.05 
At* 0.001 

^equilibration ^500 

f* ^0 

''thermostat ou 

''production 

TABLE II: Parameters used in the numerical simulations of this work. -/V a t ms is the number of 
atoms, Zg is the cutoff in the potential, At* is the molecular dynamics time step, ^equilibration * s 
the equilibration time, ^thermostat ^ s the time between velocity rescalings during equilibration, and 
^production * s the total length of the molecular dynamics trajectory. We have also performed runs 
with iV atoms = 32 and 256. 

V. RESULTS 

A. Equation of state 

In this section we present our results for the equation of state. We also make comparisons 
to the available experimental data and to analytical results in the high temperature limit. 
The pressure in the limit nX 3 <C 1 is given by the virial expansion 



The second virial coefficient is wel 



1 - b 2 (^f) + (Abl - 26 3 ) + 0((nX 3 ) 3 ) . (28) 

Z>2 = 3/(4a/2), and the third virial 



■known 



coefficient has been computed in 26r428j]. 63 = 0.29095295. 

In Fig. [T] we show the pressure normalized to nT as a function of the dimensionless density 
nX 3 and the temperature in units of the Fermi temperature, T/Tp. The data points show 
the results of the molecular dynamics simulations. The results are compared to the virial 
expansion at second and third order (dotted and dashed lines), and to a parameterization 2 
(band) of the experimental data from [30|. More accurate results for the equation of state 



have recently been published by the MIT group 311] • In the range of temperatures that are 



See appendix A of 29(. 
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FIG. 1: Pressure in units of the ideal gas pressure as a function of T/Tp (bottom axis) and nA 3 (top 
axis). The data were obtained from quasi-classical molecular dynamics simulations with different 
numbers of particles in a periodic box. The dashed and dotted lines show the virial expansion at 
second and third order, and the band is a parameterization of the experimental data of the ENS 
group [3C|. 



of interest to us, T/Tp > 0.5, the more recent data agrees with the earlier data within the 
errors indicated by the thickness of the band. 

We observe that the data follow the second order virial expansion for nX 3 < 0.2. This 
agreement is of course a consequence of the way the potential was constructed, but it serves 
as a useful check of the molecular dynamics (MD) simulation. We also observe that the full 
MD simulation is better behaved than the virial expansion. While the virial expansion is 
not useful for n\ 3 > 0.3 the MD results follow the data up to densities nX 3 ~ (0.5 — 1.0). 
The MD results do not reproduce the rapid increase in the pressure for nA 3 > 2. Based on 
the sign of the third virial coefficient it is reasonable to assume that the MD results in this 
regime could be improved by including a repulsive three body potential. 



10 



0.2 0.4 0.6 0.8 1 

rX 



0.2 0.4 0.6 0.8 1 

rX 



FIG. 2: Radial distribution functions and extracted from molecular dynamics simulations 
with 108 particles performed at different values of the diluteness parameter nA 3 . The solid lines 
show the analytic result in the high temperature (low density) limit. The right panel shows the pair 
correlation function for unlike spins scaled by r 2 . The intercept is proportional to Tan's contact 
density C. Note that the correlation function for rX < 0.05 is sensitive to the regulator in the 
potential. The thin solid lines show the fit described in the text. 



B. Pair correlation function 



For a homogeneous system the pair correlation function is denned as 

V 

N(N 



G(r, t) = N( ^_ 1) |r1(0) - fj(t)\ ] ) . (29) 



The pair correlation function measures the probability of finding two particles separated 
by a distance r and time t. For t = the quantity G(r, 0) is also known as the radial 
distribution function, or as the Fourier transform of the static structure factor. For a two 
component system we can define two correlation functions, G^ = and = G^, which 
characterize the probability of finding two particles of the same or opposite spin close to one 
another. 

Fig. [2] shows the radial distribution function for a system of 108 particles (A^ = iVj. = 54) 
at different densities. The solid curves show the high temperature (n\ 3 1) limit 

u eff (r) 



G(r, 0) = exp 



(30) 



k B T 

where u e ff is given in Eq. ( 124")) . We observe that the like spin correlation function is repulsive, 
which is a reflection of the Pauli principle, and the unlike spin correlation function is strongly 
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FIG. 3: Contact density C in units of NkF as a function of nA 3 and T/Tp. We also show the 
high temperature limit from the virial expansion at second and third order, and the result of two 
T-matrix approximations, labeled GPF/GqGo (see text). 



attractive, which is a consequence of the attractive interaction in the spin singlet channel. 
The range of the correlation function is equal to the thermal wave length. The MD results 
show that at non-asymptotic temperatures the correlation function are more short range. 
We also observe that the equal spin correlation becomes attractive at intermediate range, 
and that there is less attraction in the opposite spin channel. The pair correlation function 
at zero temperature was studied using Green Function Monte Carlo (GFMC) 0, 32)]. At 
T = the pair correlation function has the same shape as in the high temperature limit, 
but the range is set by hp 1 . 

The right panel of Fig. [2] shows the short distance behavior of the unlike spin correlation 
function. Tan observed that the correlation function is proportional to 1/r 2 at all tempera- 



tures 



33 



34|. The 1/r 2 term is governed by a universal parameter known as Tan's contact 



density C, 



C 



1 



2 

ar 



(31) 



This relation combined with 
high temperature limit 35^37 1 



]q. ( 130]) shows that the contact density scales as T 1 in the 



C 



mT 



(32) 



12 




FIG. 4: Unequal time pair correlation function for nA 3 = 1.0 (left) and nA 3 = 3.5 (right) at At = 
0.2, 0.4, 0.6, where At is expressed in units of Xy/m/T. The solid curves show the corresponding 
At = result from Fig. [2j 



Fig. |2] shows that at non-asymptotic temperatures the contact is smaller than the limiting 
form given in Eq. ( 13"2"j) . The pair correlation function in the limit r — > is sensitive to 
the regulator Iq in the potential. In order to extract the contact density we fit the unlike- 
spin correlation functions at short distances (Iq < rX < 0.5) with the functional form 
(r\) 2 G^i(r\) = C/(4:Tr 2 n 2 ) + a (rX) p , where a ,p and C are treated as fit parameters. The 
quality of the fit can be seen from the thin solid lines in the right plot of Fig. El The value 
of the contact density obtained from the fit is shown in Fig. The size of the data points 
approximate the error in the fitted contact density estimated by varying the size of the fit 
region by 10% in either direction. 

In Fig. [3] we also compare the MD result for the contact density to recent T-matrix 
calculations. The line labeled GPF/G G shows the results for two particular truncations, 



the Generalized Pair Fluctuation (GPF) theory of Nozieres and Schmitt-Rink 37|, |38(, and 



the non-self consistent (GoGo) T-matrix approximation of Palestini et al. 3 41] . We find very 
good agreement between the results of the MD simulation and the T-matrix calculations for 
all densities we have studied, 0.05 < nX 3 < 2. We also note that the T-matrix calculations 



We have taken these results from the compilation in |37|. The GPF and GqGq result s ag ree within the 



width of the band in Fig. [3] The self consistent T-matrix calculation of Punk et al. 39( is about 10% 
higher at T/Tp = 1. The recent Path Integral Monte Carlo calculation of Drut et al. 01 indicates that 
C/^NUf) reaches a maximum of ~ 3.4 at T/Tp ~ 0.4 and then decreases slightly to ~ 3.0 at T = 0. 
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(convoluted with suitable density profiles for finite traps) were shown to agree with the 



recent data reported in [421 ]. see [37] 



Fig. H] shows the pair distribution function at non-zero time difference. This is the first 
quantity in this work that is not directly amenable to quantum Monte Carlo studies. The 
Fourier transform of G(r, t ) yields the dynamic structure function, which has been studied 

nn 

extensively in a variety of many-body theories [1|, [2J. We have not attempted to perform 
the Fourier transform, since this would require high statistics data on a very fine mesh. 
We observe that the typical correlation time in our data is of order one in simulation units 



Xy/m/T. 



VI. SHEAR VISCOSITY 



We compute the shear viscosity coefficient using the Green-Kubo relation 

v r°° 

r)(xy) = ^f] o ( /, ^ / )/ , ,, / (0)) dt , (33) 
where P xy (f) is the stress energy tensor evaluated at time t |43l. l44j| . 



N N N 



p xy (t) = m + \ ( Xi ~ x ^ ( yi ~ y ^ ■ ^ 

i=l i=l j^i 

There are five independent stress tensor correlation functions that can be used to determine 
the shear viscosity. We denote the corresponding estimates by r)r xy \,r)(y Z \,r]r xz ' ) ,'r]f xxy y' ) and 
V(yyzz)y an d use the average of the five measurements as our final result 4 . The results for 
N = 108 and N = 256 is shown in Fig. [5j 

We can check the MD simulation by comparing the numerical result to the expected 
behavior in the high temperature limit. The calculation of the transport cross section and 
the shear viscosity is explained in the Appendix. For a two component system with the 
classical interaction given in Eq. (|24p the high temperature limit of the shear viscosity is 

= 75 ^ (35) 
ci 8(5 + 7r)nA 3 ' y ' 

This result is shown as the high temperature limit of the solid black line label 'MD Fit' in 

Fig. [51 and we observe that the agreement with the MD results in the regime nA 3 < 0.15 is 



JL 
hn 



4 T)t yz \ and r\(x Z ) can be obtained trivially from Eq. (j33]). The diagonal (xxyy) Kubo formula is given by 
V(xxyy) = t^t J™ (i p ^{t) - Py y (t)) (P xx (0) - P y y(0))) dt, and V(yyzz) is defined analogously. 
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FIG. 5: Molecular dynamics results for the shear viscosity r/ in units of HX~ 3 as a function of nA 3 
and T/Tp. The solid black line shows the fit to the MD results as described in the text, and the 
dashed blue line shows the quantum mechanical high temperature result. The band shows the 



experimental data published in 



46j |. Note that the data represent trap averages, and that we 



have used the ideal gas temperature at the center of the trap. 



very good. The classical result is about 20% larger than the (almost exact) quantum result 



in the high temperature limit 
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hn 



4l5| 

457r 3/2 



T 



3/2 



(36) 



64v/2 \T F J 16 (nA 3 ) 
This result is shown as the dashed blue line in Fig. [5j The discrepancy between the classical 
and quantum calculation is due to a combination of two effects, discussed in more detail in 
Appendix |Al The first is the fact that the unlike spin potential does not exactly reproduce the 
quantum mechanical transport cross section. The second effect is that the like spin potential, 
related to Pauli repulsion, leads to a finite transport cross section, even though there is no 
scattering in a quantum system of like spins interacting by a pure s-wave potential. 



Fig. E]also shows a parameterization of the viscosity measured in |45|, |46|. The thickness 
of the band approximates the statistical errors in the measurement. Note that the measured 
viscosity is a trap averaged quantity, and as a result there is a ~ 20% discrepancy between 
the data and the theoretical result in the high temperature limit (dashed line). The very close 
agreement between the measured viscosity and the MD simulations is therefore somewhat 
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of a coincidence. 

We note that the molecular dynamics simulation reproduces the A -3 scaling of the shear 
viscosity. We also emphasize that the numerical discrepancy between the classical and 
quantum result is quite small. It is therefore reasonable to extract an estimate of the 
leading density dependence of the shear viscosity from the MD simulation. We have fit the 
MD results shown in Fig. [5] with the expression 

^ = ^(l + c 2 (nA 3 )) , (37) 

where rjo was fixed using Eq. (1351) . We find C2 — 0.32. The fit well represents the MD results 
in the density range studied in this work, 0.1 < (n\) 3 < 2, including the rise in 1]/X 3 around 
(nA) 3 ~ (0.5 — 2.0), which is also seen in the experimental data. 



VII. DISCUSSION AND FUTURE WORK 



We have presented a new method for studying the dynamics of cold atomic gases based 
on molecular dynamics simulations with an effective quantum potential. In this work we 
have restricted ourselves to two-body interactions. In this case the MD simulation exactly 
reproduces the second virial coefficient and the pair correlation function in the dilute limit. 
We have measured the equation of state, the pair correlation function, and the shear viscosity 
for a range of densities 0.1 < nX 3 < 2.0. We find that we can reproduce the experimentally 
measured equation of state for densities nA 3 < (0.5 — 1.0). This is an improvement over the 
virial expansion, which breaks down nA 3 > 0.3, and we expect that the range of applicability 
can be extended by including a three-body force. 

We have also measured the static and dynamic pair correlation functions as well as the 
shear viscosity. We find that in the dilute limit the contact C scales as ra^n^/(mT), in 
agreement with the prediction in 3^, [3? | . Higher order correlations suppress the contact 
relative to the asymptotic behavior. Excellent agreement between the MD simulations and 
T-matrix calculations is seen for nA 3 < 2. 

In this regime the temperature dependence of the shear viscosity is well described by the 
functional form rj = r/oA~ 3 (l + C2nX 3 ). The scaling of rj with A~ 3 agrees with the expected 
behavior at unitarity. The parameter 770 differs from the exact quantum mechanical behavior 
at high temperature by a factor 1.2. We discuss the origin of this factor in the appendix, 
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but it would clearly be desirable to understand the physical origin of the discrepancy more 
clearly, and to investigate whether there are any improvements of the quasiclassical MD 
method that reproduce the correct transport cross section. The main result of the MD 
simulation is that the density dependence of the shear viscosity is weak, — 0.32, and that 
it tends to increase the shear viscosity. 

There are a number of additional applications or further extensions of the work presented 
here: 

1. BEC-BCS crossover: In this work we focused on a dilute Fermi gas at unitarity. 
Clearly, the same methods can be used to determine the effective potential and the 
leading nX 3 correction to pair correlation functions and transport coefficients for the 
full BEC-BCS crossover. 

2. Three-body forces: The method discussed in this work can be systematically improved 
by including ra-body forces (n > 3). In the case of three-body forces this appears 
quite tractable. The three-particle Slater sum can be computed in hyper-spherical 
coordinates, and molecular dynamics with three-body interactions is computationally 
feasible. 

3. Other transport coefficients: The methods used in this work can also be used to 
measure other transport properties like the spin diffusion constant studied theoreti cally 



in [47j and recently measured by the MIT group [48j], or the thermal conductivity [49 ]. 
It would be interesting to determined whether nX 3 corrections are universal. It is also 
interesting to determine the bulk viscosity as one goes away from the unitary limit 
a —7- 00. At unitarity the bulk viscosity is expected to vanish based on theoretical 
reasons 



50] . and existing experiments are consistent with these arguments 



511. 



4. Lower dimensional systems: Recently Vogt et al. measured the shear viscosity of a two- 



dimensional Fermi gas 52[. In two dimensions fluctuations are expected to be more 
important than in three dimensional systems. These effects are difficult to include in 
kinetic theory, but are automatically included in molecular dynamics. 
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Appendix A: Shear viscosity of dilute classical and quantum gases 

In this appendix we summarize the calculation of the shear viscosity for a dilute classical 
gas interacting via the potentials given in Eq. ( 1241) . This result provides an important check 
for the calculation of the viscosity using the Green-Kubo formula. For both quantum and 
classical gases the calculation of the shear viscosity using the Chapman- Enskog method jo^ ] 
yields 



5 yJirnik B T 

" = 8^^< (A1) 



where a tr is the energy averaged transport cross section 

poo 

a tI = / 7 V^ ( r tr (7) d 7 , (A2) 
Jo 

and Otr (7) is the transport cross section. The parameter 7 = y 1 Ejk B T is the dimensionless 
relative velocity of the two-particle system. The calculation of the transport cross section is 
different in the classical and quantum case, both because the cross section is represented in 
a different way, and because of the presence of symmetry factors in the quantum mechanical 
calculation. 

1. Dilute classical gas 

In the classical case the cross section is computed from classical trajectories in the po- 
tential. We can write cr tr (7) as a one dimensional integral over the impact parameter b 



C"tr (7) 



poo 

2tt / [1 - cos 2 x(b,E)] bdb. (A3) 
Jo 

In this expression the scattering angle x{b, E) is a function of the impact parameter and of 
the kinetic energy in the center-of-mass system through the relation 

* ( ^ = '- a iL i-^-gM/* ' (M) 
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where U(r) is the interaction potential and ro is the distance of closest approach. The 
parameter r$ is determined by the solution of 

1 - b 2 /r 2 Q - U(r )/E = . (A5) 

For the potentials of interest, Eq. (|24j) . all of the above expressions must be evaluated 
numerically. We find the following results for a two particles interacting through the 
and Uft potentials, 

<?£ = h 2 x (1.004) , (A6) 

5? = T^A 2 x (1.0006) , (A7) 
15 

where the factors 1.004 and 1.006 are the results of numerical integrals. This work studies 
a two-component gas with equal numbers of spin up and down particles. In the dilute limit 
this system can be treated as a classical mixture and the effective transport cross section is 
the average of the a t f and cross sections. The shear viscosity is 



'tr allKJ - w tr 

7] 75V2tt 
hn 



ci 8(5 + tt) n\ 3 ' 

which agrees with the high temperature MD simulation as shown in Fig. [5j 



(A8) 



2. Dilute quantum gas 

The quantum mechanical cross section can be represented in terms of the scattering 
phase shifts. In quantum mechanics we also have to take into account the symmetry of 
the wave function, which depends on whether the particles are distinguishable or not. For 
distinguishable particles the transport cross section is 

(t)^E {l \21f3) 2) sin2 [6l+2{k) ~ m] ' (A9) 

where k = \JmE/h with E the center of mass energy. The parameter 7 is given by 7 = 
a/ E/k B T as above. For indistinguishable particles 

(7) = f E {l+ ( 2l + 3) 2) Sin2 [6l+2(k) - m] ' (A10) 
Z=e/o 

where the sum is restricted to even/odd (e/o) angular momenta for bosons/fermions. For a 
two-component Fermi gas interacting via an s-wave interaction we can treat the collisions 
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between atoms of opposite spin as occurring between distinguishable particles. At unitarity 
the transport cross section is 

<#(7) = |£, (AH) 
and the energy averaged cross section is given by 

= ^ 2 . (A12) 

We find that the quantum cross-section is a factor of two larger than the analogous classical 
cross-section, Eq. (IA6[) . As in the classical case, the effective transport cross-section entering 
into the viscosity consists of an average over the anc i a tt channels. As the latter is zero 
the overall cross-section is reduced by a factor of two. The final result is 



hn 16(nA 3 ) 



(A13) 
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